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Abstract: We propose a method to reduce the computational effort to solve a partial differential 
equation on a given domain. The main idea is to split the domain of interest in two subdomains, 
and to use different approximation methods in each of the two subdomains. In particular, in one 
subdomain we discretize the governing equations by a canonical scheme, whereas in the other one 
we solve a reduced order model of the original problem. Different approaches to couple the low- 
order model to the usual discretization are presented. The effectiveness of these approaches is tested 
on numerical examples pertinent to non-linear model problems including the Laplace equation with 
non- linear boundary conditions and the compressible Euler equations. 
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Methodes Iteratives pour la Reduction de Modeles par 
Decomposition de Domaine 

Resume : On propose une methode pour reduire les efforts de calcul pour resoudre une equation 
aux derivees partielles sur un domaine donne. L'idee principale est de diviser le domaine considere 
en deux sous-domaines, et d' employer differentes methodes d' approximation dans chacun des deux 
sous-domaines. En particulier, dans un des sous-domaines I'equation en question est discretisee par 
une methode canonique, tandis que dans F autre un modele d'ordre reduit du probleme original est 
utilise. Des strategies differentes pour coupler le modele d'ordre reduit a la discretisation habituelle 
sont presentes. L'efficacite de ces approches est testee sur des exemples numeriques pertinentes pour 
des problemes modeles non lineaires, notamment I'equation de Laplace, avec des conditions limites 
non lineaires, et les equations d'Euler compressibles. 

Mots-cles : modeles reduits, decomposition de domaine, ecoulements compressibles 
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1 Introduction 

In this contribution we are concerned with the coupUng between a full order simulation and a reduced 
order model. The idea is to reduce the extent of the domain where we perform a canonical numerical 
simulation by introducing a low-order model which describes the solution far from the region of 
interest. By reducing the extent of the domain we aim at reducing the costs in terms of required 
memory as well as in terms of computational time. 

In a broad sense, there exist many applications where far from the boundary the solution is 
weakly dependent on the details of the boundary geometry. In such regions we use a reduced order 
model based on proper orthogonal decomposition (POD) [:4| to solve the problem. This approach 
allows a representation of the solution by a small number of unknowns that are the coefficients of 
an appropriate Galerkin expansion. Therefore away from a narrow region close to the boundary of 
interest the number of unknowns to be solved for is drastically reduced. This idea was previously 
explored in the context of transonic flows with shocks |]3], El- Here we extend those works by 
adapting to that context some classical domain decomposition techniques. 

In the following we discuss three possible methods to do that. The first is based on a Schur 
iteration where the solution of the low-order model is obtained by a projection step in the space 
spanned by the POD modes. The second is in the same spirit but instead of a Dirichlet-Neumann 
iteration we employ a Dirchlet-Dirchlet iteration in the frame of a classical Schwarz method. The 
last approach is of different nature since the solution of the low-order model is not simply based on 
a projection in the space of the POD modes. It takes into account in a weak sense the governing 
equations by minimizing the residual norm of the canonical approximation in the space spanned by 
the POD modes. 

The numerical demonstrations shown in the following are relative to two models: the Laplace 
equation with non-linear boundary conditions modeling radiative heat transfer and the compressible 
Euler equations in a nozzle. Since this method can be of interest for optimal design applications, 
where many different geometries must be tested to improve performance, in some cases we have 
explored the idea of simulating by usual discretization methods just the region where the geometry 
changes, modeling the rest by POD. 

Like all other approaches based on POD, a solution database is necessary to build the basis 
functions, therefore this method will be useful when many computations for relatively similar cases 
are to be performed, like for example in shape optimization, see 0] . 

2 Solution by projection of the solution trace in the space spanned 
by the POD modes 

2.1 Approximation of the Steklov-Poincare operator by POD 

In order to explain the method we take a particular case. Let us consider the Laplace equation Au = 
definedjnside a square n = [0, 1] x [0, 1]. Let d € [0, 1], r^i = [0, d] x [0, 1], = [d, 1] x [0, 1] 
and r = ill nri2 die interface between the two sub-domains, see fig.[I] We have Dirichlet conditions 
on the right boundary (ufl) as well as on the upper (uu) and lower (uu) boundaries. 
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r= QinQ2 



1 -d 



Figure 1: Problem set-up: Schur 



We want to solve this problem for different values of the Dirichlet data on the left (ul) boundary. 
To that end we build an appropriate solution database for a given set of boundary conditions on the 
left side. In particular, the Dirichlet data on the left boundary is denoted by g^, 1 < k < N. 

Let the functions be the discrete solutions of the boundary value problem posed in CI 
for different k, and let U2 be an harmonic function, restricted to 1^2, defined as follows: U2 = 
^/^J2k=i '^n^ ■ compute a Galerkin base of the form (pi — J2k=i ^ikiu^^ ~ U2) where u"^^ 
is the restriction of to Vt2- The coefficients 5;^ are found by POD as explained in 18]. It can be 
shown that this base gives by construction an optimal representation of the original data set . 

Let us define M2 = U2 + X^fii (^i4>i^ where M is much smaller than the number of discretization 
points in Vl2- For an arbitrary Dirichlet condition on the left boundary of 51, we want to determine 
the discrete solution by a canonical approximation in 51i and by the above defined Galerkin repre- 
sentation in Vl2- 

One simple way to implement this idea is to solve the problem by Dirichlet-Neumann iterations. 
To this end, we follow the steps below: 

1. solve the problem in ili by any discretization method (FD, FEM, etc.), imposing Neu- 
mann b.c. on F; 

2. on interface F project the trace of the above solution in the subspace spanned by the 
traces of the POD modes 0i; 

3. recover U2 as the prolongation of the trace of U2 on F inside O2 by using the POD 
modes; 
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4. set dui/dn = du2/dn on T\ 

5. goto (1) until convergence is attained. 

This is just one possible solution algorithm, corresponding to a classical domain decomposition 
method (Schur complement). Another approach consists in solving the problem all at once, as 
detailed in the following. Let us define Ai the discretized operator acting on ui, the restriction of 
the unknowns belonging to fii; Ay the discretized operator acting on wr, the unknowns belonging 
to r and A2 the discretized operator acting on U2, the restriction of the unknowns belonging to Vl2- 

The discretized non-linear problem in il can be written 

Ai Bi 

B{ Ar I I = I /r I (1) 

B2 A2 

where Bi and B2 are appropriate interface matrices and /i, fr, f2 take into account the boundary 
conditions. 

From ([T]i we have 

Aim + Biur = fi 
B{ui + (Ar - BlA^^B2)tiT = fr - B\A^^ f2 (2) 

The matrix Ay — i?2^^^i32 is the discrete counterpart of the Steklov-Poincare operator for ^2, see 
||6l . Consider now the second step of the solution algorithm proposed above. Let a G be a 
vector of components oi . . . om and c e R^^ a vector of components ci . . . cm- Posing the trace 
of (t>i on r, we take 

M 

Mr - Ck^Pk 



arg mm 



k=l 



(3) 



where ||-|| is the norm induced by the canonicaH^ scalar product, noted by (•,•)■ 

Solution of ([3]) reduces to the solution of the linear problem X]i=i i'Pij Vj) = ("r, 'Pj), 1 < 
j < M. Therefore ~ {ur, Pi), where Pi = 12jLi [{fi: fj is a constant vector computed 

once for all from the POD modes. 



At this point we approximate U2 with U2 and substitute in ([TJ. Since B\u2 = -S2U2+X]f=i Oi-S* ' 



we have B2U2 = B\u2 + X]i=i B24>i{uT, Pi)- Finally, letting S2 = X]f=i B2<j)iPi we obtain the 
approximation of (|2]i 

Blm + [At - S2)ur = fr - Blu2 (4) 

where 52^2 = ^2^2^^/2- Matrix 5*2 is the approximation of the discrete Steklov-Poincare operator 
obtained by the POD expansion. Equations (|2|i can of course be solved simultaneously by a standard 
linear solver 

Just like for the usual Steklov-Poincare operator, (|4| amounts to a non-local boundary condition 
for the problem posed in ili . The main advantage of this approach compared to computing explicitly 
5*2 is that we do not need A^ ^ to build S2 ■ 
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In the following we present some numerical applications. The first case is a plain application of 
the Schur complement as it was explained above to a non-linear case. The second case is based on 
the Schwarz method. The third numerical experience is a variant of the all at one method. 

2.1.1 Schur complement 

A second order finite differences (FD) method coupled to a fix point iteration is used to solve the 
Laplace equation inside the square domain shown in fig.[T] The left Dirichlet boundary condition is 
varied to build the needed database. In particular, ul — sin {kny) + y, I < k < 49. The boundary 
conditions on the other sides are, referred to fig.[T] the following: on uu : u = 1, on ud '■ u = and 
onuR:u^-4 + §^^^0. 

The domain is split at d = 1/3. Then, the POD basis functions are generated on using the 
previously computed database. In order to check the accuracy of the method, a boundary condition 
which was not included in the database used to build the POD modes is imposed on the left boundary: 
Ul = y^, and a second order FD method is used to solve the problem in fli. We use 6 POD modes 
to recover U2 inside il2- Figure |2] presents the result of the test by means of the distribution of the 
relative error between the FD solution on the entire domain and the approximate Schur complement 
approach. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 O.B 0.9 I 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 I 



(a) (b) 

Figure 2: Distribution of the relative error between the solution obtained by the present method and 
the solution obtained by a second order FD methodon the whole domain, fli (a), fl2 (b). 



2.1.2 Schwarz method 

In the following, a convergent-divergent domain il is considered. As before, is divided in two 
subdomains, 51i and ^2, in a way that there exist an overlap region ftov = fii n shared by both 
subdomains (see fig.O. The Laplace equation is solved with the boundary conditions detailed below 
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to generate a database of k solutions with 1 < fc < 60, this time varying the geometry of The 
solution is obtained by means of the finite element method as implemented in [l5l, using PI elements 
on a triangular non-structured mesh and a fixed point iteration. The boundary conditions imposed 
are: on it^ : u = i • on ujy : u ~ 1, on u/j : m = and on ur : — + ^ = Q. 




Figure 3; Problem set-up: Schwarz 

For a geometry of the divergent part which is not included in the database, the solution is deter- 
mined following a similar approach to that described for the Schur complement but employing this 
time the classical Schwarz method (see, for example, 16]): 

1 . solve the problem in Oi by any discretization method imposing Dirichlet b.c. on r,\; 

2. on project the trace of the above solution in the subspace spanned by the trace of 
the POD modes 4>i; 

3. recover U2 as the prolongation of the trace of U2 on r^^ inside O2 by using the POD 
modes; 

4. set ui = U2 on r,\; 

5. goto (1), n = n + 1, until convergence is attained. 

Four POD modes are used to recover U2 inside Figure |4] shows the results obtained for this 
case, again in terms of the relative error between a numerical solution (FEM PI) on the entire domain 
and the approximate Schwarz method. 
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2.1.3 Newton method 

In this case the two-dimensional compressible Euler equations are solved 
= —V • grad S 

= —V- grad a — ^2^a div V (^^■^ 
= - grad V ■ V - a (^::^ grad a - a grad 5^ , 

with 7 = 1.4 for air and being S the entropy, a the speed of sound and v = [u,v] the velocity, in a 
convergent-divergent nozzle on a structured mesh. We consider shockless flows. 

The A-scheme ||9l is used to solve the equations. Total temperature, total pressure and the flow an- 
gle are imposed at the inlet; static pressure at the exit and impermeability at the walls. The complete 
system is solved by Newton iterations. The resulting linear problems are solved by preconditioned 
GMRES iterations |7J. 

Using this code, a database of 90 snapshots is computed. The geometry of the divergent part 
of the nozzle is changed, fii in fig. [3] Furthermore, the static pressure at the exist is also varied 
taking uniformly spaced values in the interval p = [0.94, 0.99] with step 0.01. This corresponds 
to 15 snapshots for each pressure step. Since the flow is shockless, total temperature as well as 
total pressure are constant across the nozzle. Therefore, in order to solve in fii, the only unknown 
boundary condition is the flow angle on F^. Hence, a low-order representation of the flow angle v /u 
for the convergent part, i.e., ^l2, is constructed retaining 10 POD modes. 
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We consider a case for which the geometry of the divergent part of the nozzle and the static 
pressure imposed at the exit do not belong to the set used to build POD database. The two sub- 
problems are coupled by a least squares approximation on the overlapping domain. This is just a 
variant of the Schwarz method presented in the previous section. The only difference is that the two 
coupled problems are solved all at once by a Newton method. 

The relative errors restricted to ili are reported in table[T]for some of the flow variables. These 
errors are computed with respect to the numerical simulation on the entire domain, and are quantified 
in terms of relative error in norm, i.e., the norm of the difference between the solution obtained 
by the present method and the reference solution divided by the norm of the reference solution. 



Variable 


e{u)% 


e{v)% 


e{Ma)% 


e{v/u)% 


Error 


1.1 


2.25 


1.14 


0.75 



Table 1: Relative percentage errors (in norm) for the flow variables: u (horizontal velocity), v 
(vertical velocity). Ma (Mach number) and v/u (flow angle). 

Figure |5j a) presents the distribution of the Mach number obtained by the full numerical simula- 
tion on the entire domain. The same quantity obtained when applying the proposed method is shown 
in fig. |5jb). Finally, fig. |5jc) shows the distribution of the relative percentage error. In figure |6]the 
flow angle v /u is considered. 

3 Solution by minimization of the residual norm in the space 
spanned by the POD modes 

An alternative way to couple the low-order model to a detailed simulation is to look for an approx- 
imate solution in the reduced order function space that takes into account the governing equations. 
Hence, the main difference with respect to the approach described in section |2] is that the approx- 
imate solution is found by minimizing the residual norm of a given discretization scheme on the 
whole domain O2 rather than by projecting the trace of the solution in the space spanned by POD 
modes. 

Again the total pressure and the total temperature are constant across the nozzle and therefore we 
can write any other variable as a function of the local Mach number and the ratio v /u. In particular 
let U be the array of the couples {Ma, v/u) for all the grid points belonging to il2- We start by 
representing this vector in the original Q-dimensional discrete space by a linear combination of 
basis functions, U = '^f^i oa^i, with M <C Q- The arrays <I>i have the same structure of U and 
they have been obtained by POD. The idea is to satisfy the compressible Euler equations (|5]l in a least 
squares sense over ^2- In other words, let E{a) be the discrete residual of the governing equations 
as a function of a = (ai , . . . , om) and let 

lia) = ^E^{a)Eia) (6) 
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(a) 




(b) 




(c) 

Figure 5: Distribution of Mach number: (a) Numerical solution (X-scheme). (b) Present method, 
(c) Relative percentage error distribution on fli. 



be the residual norm. The solution in the POD function space a* is found by setting 

a* — arg minQ,/(a) 



(7) 
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1 1.2 1.4 1.6 



(a) 




(b) 




(c) 

Figure 6; Distribution ofv/u. (a) Numerical solution ( X-schetne) on the entire domain, (b) Present 
method, (c) Relative percentage error distribution on fli. 



This is equivalent to a system of non-linear equations 

dl _ dE^ia*) 



dai 



E{a*) = (a*)E{a*) = Vi £ [1, . . . , M] 



(8) 
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where is the system's Jacobian. 

Equation dH) is solved by a quasi-Newton method where we take 

E{a*) = E{a°) + J{a°)Aa + 0{Aa'^) (9) 

and after substituting in (O we obtain 

J^JAa, = -J^E (10) 

Therefore, we are left with the solution of equation dTol ) at each step of the quasi-Newton algorithm 
employed to solve (O. 

Using the same database generated in Section 12.1.31 two low-order basis restricted to the con- 
vergent part of the nozzle are constructed. One for the flow angle v ju and the other for the Mach 
number. Only 5 POD modes are retained for both the flow angle and the Mach number 

This method is tested for the same case as for the Newton method presented in Section 12.1.31 
The only difference is that the solution in ^2 is found by minimization of the residuals norm in the 
space spanned by the POD modes. Otherwise, the full numerical simulation is employed on fii, and 
a classical Schwarz overlapping method is used to iterate the solution to convergence. 

The reference Mach number on the entire domain, the Mach number obtained by the present 
method on fii as well as the distribution of the relative percentage error on fJi are shown in figures 
[8l;a),[8l:b) and[8tc), respectively. 

In fig. |7] the results for the flow angle vju are shown while table |2] presents the relative errors 
obtained in ili. 



Variable 


e(u)% 


e{v)% 


e(Ma)% 


e{vlu)% 


Error 


0.05 


0.26 


0.05 


0.13 



Table 2: Relative percentage errors (in 1? norm) for the flow variables: u (horizontal velocity), v 
(vertical velocity). Ma (Mach number) and vju (flow angle). 

These results show a much better accuracy as compared to those of section 12.1.31 as one can 
conclude by comparing tables[T]and|2] 

Finally, in order to asses sensitivity of the results with respect to the minimization problem, we 
show in table[3]the error in ^1 when we use for a the initial guess. In table |4]we present the error 
at the end of the Schwarz iteration. It is seen that the fact of approximating the solution of the Euler 



Variable 


e(u)% 


e(v)% 


e(Afa)% 


e{vlu)% 


Error 


35.52 


45.65 


39.28 


6.25 



Table 3: Relative percentage errors (in 1? norm) for the flow variables: u (horizontal velocity), v 
(vertical velocity). Ma (Mach number) and vju (flow angle). 

equations with appropriate boundary conditions, even in a crude way, in significantly improves 
the solution with respect to the initial guess. 
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(a) 



(b) 




(c) 



Figure 7: Distribution ofv/u. (a) Numerical solution ( X-schetne) on the entire domain, (b) Present 
method, (c) Relative percentage error distribution on fli. 
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(a) 



(b) 




(c) 

Figure 8: Distribution of Macxh number on (a) Numerical solution (X-scheme). (b) Present 
method, (c) Relative percentage error distribution on fli. 

4 Discussion 



The fact of reducing the extent of the computational domain does not guarantee that one can get a 
solution faster as compared to solving the problem on the entire domain. In particular, let us consider 
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Variable 


e(ii)% 


e(«)% 


e(Ma)% 


e{v/u)% 


Error 


2.01 


3.94 


2.08 


2.90 



Table 4: Relative percentage errors (in norm) for the flow variables: u (horizontal velocity), v 
(vertical velocity). Ma (Mach number) and v/u (flow angle). 

the method described in section 12.1.31 It would in principle be the less demanding in terms of 
computational time since it just amounts to a non-local boundary condition in the frame of a Newton 
method. However, this is not necessarily the case since the convergence of the implicit iteration is 
spoiled by such non-local boundary condition and the number of Newton steps to attain convergence 
goes from 7 on the whole domain, to 1 1 when we only solve on Oi. The slower convergence rate 
may lead to comparable costs in terms of CPU time, as a function of the fraction of the domain that 
is actually resolved. On the other hand, the fact of reducing the number of grid points is reflected 
almost proportionally on the memory requirements. However, the Jacobian matrix will have a non- 
sparse block corresponding to the non-local boundary condition induced by the approximation of 
the Steklov-Poincare operator 

Concerning the method described in section[3]the cost of each step of the proposed minimization 
algorithm can be split up as follows, i) The computation of J. The Jacobian is evaluated by one- 
sided finite differences. Consequently its cost in terms of CPU time and memory requirements is 
proportional to the dimension of the low-order space M and the cost of computing the residual 
vector, i.e., cost(J) = M x cost(£'). Since Af is 0(10), the Jacobian matrix can be formed after 
few residual evaluations that are cheap in terms of CPU time. Memory requirements could become 
prohibitive for very big problems only on serial architectures, ii) The computation of the symmetric 
matrix A = J^J. The dimensions of A are M x M, with M = 0(10). The required floating point 
operations are {M x -A/)/2 x N and they are less than the operations needed for the computation 
of E. iii) The cost of finding the solution to the linear system of size M x M can be neglected. 
Summarizing, it can be deduced that the cost is equivalent to some residual evaluations. The number 
of iterations needed in order to find the minimum are generally less than five, yielding an overall 
cost negligible with respect to canonical CFD calculation on il2- 

We remark two major limits of these approaches. The first is of course that the results depend 
to a large extent on the database used for the POD modes. If the configuration under consideration 
lies in a region of the parameter space far from that explored when building the database, then the 
approximation error can be large. It is true, however, that the low-order model should be used 
where the solution does not strongly depend on the boundary conditions or the geometry. Another 
limitation is that an efficient way to improve the approximation quality comparable for example to 
grid refinement is not available. In principle, we would like to increase the approximation accuracy 
by enriching the functional space in which the solution is sought, based on some objective criteria. 
Unfortunately a general framework for such improvement is not presently available and it is the 
object of present research. 

In conclusion we presented some possible implementations of a method to reduce the extent 
of the computational domain in the numerical solution of partial differential equations. The idea 
of using models that take into account different physical phenomena in different subdomains is 
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old. Here we revisited this approach using a low-order model in the framework of classical domain 
decomposition techniques. The results in terms of the approximation error are promising for all of 
the cases that we showed. A major challenge to be pursued is to find a viable a posteriori error 
estimation technique for iteratively adapting the POD function space. 
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